global name "s2_dokrig"

global datafile "datacvr" 
global datafilemi "mult_impu"


global main "C:/Users/silvio/Documents/CVR/ryp"
global output "$main/output"  
global data "$main/data"  
global code "$main/code"  
global excel "$main/excel"  
global odata "$data/Intermuestra1-2"  

global coord "$data/coordinates"  
global geocoord "$data/geocoord"  


*1 Fix DGP. Simulation with parameters and best model from best data for SL
cd $main
capture log close
log using "$output/${name}.log", replace

matrix drop _all
set more off

set matsize 1392
scalar nstrata=58                      /*=58, CVR strata*/
scalar signi=2


use "${data}/${datafile}" , replace


sort i perpe k
ren y yy
mkmat yy

forval perpel= 1/3{            /*Original*/
scalar perpeh=`perpel'
do "$code/cvrkrigrun"
matrix list mao
matrix list sum_mao
matrix list modelmat

di "perpe:  " perpeh
di "depmean:  " depmean 
di "depvar:  " depvar
di "depvar1:  " depvar1
di "depvar2:  " depvar2
di "nlags:  " nlags
di "nstrata:  " nv
di "nestim:  " nestim
di "nestimratio:  " nestim2
di "nestime:  " nestime
di "nestimv:  " nestimv


if perpeh==1{
matrix mao_est=mao
matrix sum_est=sum_mao
matrix modelmat_e=modelmat
}
else if perpeh==2{
matrix mao_sl=mao
matrix sum_sl=sum_mao
matrix modelmat_s=modelmat
}
else if perpeh==3{
matrix mao_ot=mao
matrix sum_ot=sum_mao
matrix modelmat_o=modelmat
}
}
 /*Original*/
matrix maotot=mao_est+mao_sl+mao_ot
matrix list maotot
matrix sumtot=sum_est+sum_sl+sum_ot
matrix list sumtot
matrix sumall= sum_est\sum_sl\sum_ot\sumtot
matrix list sumall
matrix maotot=mao_est+mao_sl+mao_ot
matrix modelmat=modelmat_e,modelmat_s,modelmat_o
matrix list modelmat
matrix ones=J(1,7,1)
matrix modelmatt=ones*modelmat
matrix list modelmatt


*************************************************************************************
*	Totals
*************************************************************************************

/* Report Total */

matrix mao=sumall

matrix smao=mao[1...,3]
matrix bmao=mao[1...,2]
matrix nmao=mao[1...,1]

*matrix tmao=bmao\
matrix cil_mao=J(4,1,0)
matrix ciu_mao=J(4,1,0)
forval i= 1/4{
matrix smao[`i',1]=sqrt(smao[`i',1])
*matrix tmao[`i',1]=bmao[`i',1]/smao[`i',1]
matrix cil_mao[`i',1]=int(bmao[`i',1]-signi*smao[`i',1]+0.5)
matrix ciu_mao[`i',1]=int(bmao[`i',1]+signi*smao[`i',1]+0.5)
}


matrix mao=nmao,cil_mao,bmao,ciu_mao,bmao,smao
matrix list mao

svmat mao, names(v)
ren v1 n
ren v2 cl
ren v3 e
ren v4 cu
ren v5 em
ren v6 s

export excel n cl e cu em s using  "$excel/${name}k.xls", sheet("totals") sheetmodify firstrow(varlabels) missing(".")
clear


/* Report Strata */
foreach k in mao_est mao_sl mao_ot maotot{  /*Original*/
matrix mao=`k'

matrix smao=mao[1...,3]
matrix bmao=mao[1...,2]
matrix nmao=mao[1...,1]

*matrix tmao=bmao

matrix cil_mao=J(nstrata,1,0)
matrix ciu_mao=J(nstrata,1,0)
forval i= 1/58{
matrix smao[`i',1]=sqrt(smao[`i',1])
*matrix tmao[`i',1]=bmao[`i',1]/smao[`i',1]
matrix cil_mao[`i',1]=int(bmao[`i',1]-signi*smao[`i',1]+0.5)
matrix ciu_mao[`i',1]=int(bmao[`i',1]+signi*smao[`i',1]+0.5)
}

matrix mao=nmao,cil_mao,bmao,ciu_mao,bmao,smao
matrix list mao


svmat mao, names(v)
ren v1 n
ren v2 cl
ren v3 e
ren v4 cu
ren v5 em
ren v6 s

export excel n cl e cu em s using "$excel/${name}k.xls", sheet("s`k'") sheetmodify firstrow(varlabels) missing(".")
clear


}




/* Report regions */
run "$code/regions"
foreach k in mao_est mao_sl mao_ot maotot{   /*Origintal*/
matrix mao=`k'

matrix mao=region'*mao
scalar nregion=rowsof(mao)

matrix smao=mao[1...,3]
matrix bmao=mao[1...,2]
matrix nmao=mao[1...,1]

*matrix tmao=bmao

matrix cil_mao=J(nregion,1,0)
matrix ciu_mao=J(nregion,1,0)
forval i= 1/8{
matrix smao[`i',1]=sqrt(smao[`i',1])
*matrix tmao[`i',1]=bmao[`i',1]/smao[`i',1]
matrix cil_mao[`i',1]=int(bmao[`i',1]-signi*smao[`i',1]+0.5)
matrix ciu_mao[`i',1]=int(bmao[`i',1]+signi*smao[`i',1]+0.5)
}


matrix mao=nmao,cil_mao,bmao,ciu_mao,bmao,smao
matrix list mao


svmat mao, names(v)
ren v1 n
ren v2 cl
ren v3 e
ren v4 cu
ren v5 em
ren v6 s

export excel n cl e cu em s using "$excel/${name}k.xls", sheet("r`k'") sheetmodify firstrow(varlabels) missing(".")
clear
}


*************************************************************************************
*	Differences
*************************************************************************************

/* Report Total */
matrix mao=sumall
matrix dmao=mao[2,1..2]-mao[1,1..2]    /*d*/
matrix smao=sqrt(mao[2,3]+mao[1,3])   /*t*/
matrix tmao=dmao[1,2]/smao[1,1]    /*t*/
matrix pmao=1-normal(abs(tmao[1,1]))    /*p*/
matrix cil_mao=int(dmao[1,2]-signi*smao[1,1]+0.5)
matrix ciu_mao=int(dmao[1,2]+signi*smao[1,1]+0.5)
matrix smao=int(smao[1,1]+0.5)   /*t*/
matrix dmao=dmao,smao,tmao,pmao,cil_mao,ciu_mao
matrix list dmao

svmat dmao, names(v)
ren v1 n
ren v2 e
ren v3 s
ren v4 t
ren v5 p
ren v6 cl
ren v7 cu
format %7.4f  p
export excel n e s t p cl cu using "$excel/${name}k.xls", sheet("dmao") sheetmodify firstrow(varlabels) missing(".")
clear


/* Report Strata */
matrix dmao=mao_sl[1..58,1..2]-mao_est[1..58,1..2]    /*d*/
matrix smao=mao_sl[1..58,3]+mao_est[1..58,3]   /*v*/

matrix tmao=J(nstrata,1,0)
matrix pmao=J(nstrata,1,0)
matrix cil_mao=J(nstrata,1,0)
matrix ciu_mao=J(nstrata,1,0)
forval i= 1/58{
matrix smao[`i',1]=sqrt(smao[`i',1])
matrix tmao[`i',1]=dmao[`i',2]/smao[`i',1]
matrix pmao[`i',1]=1-normal(abs(tmao[`i',1]))    /*p*/

matrix cil_mao[`i',1]=int(dmao[`i',2]-signi*smao[`i',1]+0.5)
matrix ciu_mao[`i',1]=int(dmao[`i',2]+signi*smao[`i',1]+0.5)
matrix smao[`i',1]=int(smao[`i',1]+0.5)   /*t*/
}

matrix dmao=dmao,smao,tmao,pmao,cil_mao,ciu_mao
matrix list dmao

svmat dmao, names(v)
ren v1 n
ren v2 e
ren v3 s
ren v4 t
ren v5 p
ren v6 cl
ren v7 cu
format %7.4f  p
export excel n e s t p cl cu using "$excel/${name}k.xls", sheet("sdmao") sheetmodify firstrow(varlabels) missing(".")
clear




/* Report regions */
matrix mao_sl=region'*mao_sl
matrix mao_est=region'*mao_est

matrix dmao=mao_sl[1..8,1..2]-mao_est[1..8,1..2]    /*d*/
matrix smao=mao_sl[1..8,3]+mao_est[1..8,3]   /*v*/

matrix tmao=J(nregion,1,0)
matrix pmao=J(nregion,1,0)
matrix cil_mao=J(nregion,1,0)
matrix ciu_mao=J(nregion,1,0)
forval i= 1/8{
matrix smao[`i',1]=sqrt(smao[`i',1])
matrix tmao[`i',1]=dmao[`i',2]/smao[`i',1]
matrix pmao[`i',1]=1-normal(abs(tmao[`i',1]))    /*p*/

matrix cil_mao[`i',1]=int(dmao[`i',2]-signi*smao[`i',1]+0.5)
matrix ciu_mao[`i',1]=int(dmao[`i',2]+signi*smao[`i',1]+0.5)
matrix smao[`i',1]=int(smao[`i',1]+0.5)   /*t*/
}

matrix dmao=dmao,smao,tmao,pmao,cil_mao,ciu_mao
matrix list dmao

svmat dmao, names(v)
ren v1 n
ren v2 e
ren v3 s
ren v4 t
ren v5 p
ren v6 cl
ren v7 cu
format %7.4f  p
export excel n e s t p cl cu using "$excel/${name}k.xls", sheet("rdmao") sheetmodify firstrow(varlabels) missing(".")
clear
capture log close



** MI ****************************************************************************************************
capture log close
log using "$output/${name}mi.log", replace
scalar nstrata=58                      /*=58, CVR strata*/
scalar nmiver=20

forval perpel= 1/3{
scalar perpeh=`perpel'

matrix x2=J(nstrata, 1, 0)
matrix modelmat_m=J(7, 1, 0)


scalar diestim_mi=0
scalar nestim_mi=0

local miper=1
 while `miper'<=nmiver{
 scalar miver=`miper'

use "${data}/${datafilemi}" , replace

gen j=i
sort i perpe k

gen yy= y`miper'
*ren y yy
mkmat yy


run "$code/cvrkrigrun"
di miver
matrix list sum_mao
matrix list mao

matrix modelmat_m=modelmat_m+modelmat
scalar diestim_mi=diestim_mi+diestim
scalar nestim_mi=nestim_mi+nestim




if miver==1{
forval i = 1/58{
matrix x2[`i',1]=mao[`i',2]*mao[`i',2]
}
matrix mao_mit=mao
}
else{
matrix mao_mit=mao_mit+mao
forval i = 1/58{
matrix x2[`i',1]=x2[`i',1]+mao[`i',2]*mao[`i',2]
}
}

local miper=`miper'+1
*drop yy
} /*miper*/
 *MI****************************************************************************************************************************************************

matrix mao_mit=mao_mit/nmiver 
matrix x2=x2/nmiver 
scalar madj=(nmiver+1)/(nmiver-1)

scalar diestim_mi=diestim_mi/nmiver 
scalar nestim_mi=nestim_mi/nmiver 


scalar jhat=1
while jhat<=nstrata {
	
 	 matrix mao_mit[jhat,3]=mao_mit[jhat,3]+madj*(x2[jhat,1]-mao_mit[jhat,2]*mao_mit[jhat,2])	     /*vnctotal*/	 	 	 
	 matrix mao_mit[jhat,1]=int(mao_mit[jhat,1]+0.5)
     matrix mao_mit[jhat,2]=int(mao_mit[jhat,2]+0.5) 
	 
   scalar jhat =jhat+ 1
} /* end loop */



*MI****************************************************************************************************************************************************

matrix ones58==J(1,58,1)
matrix sum_mao_mit=ones58*mao_mit
matrix list sum_mao_mit

di sqrt(sum_mao_mit[1,3])
matrix list mao_mit
di diestim
di nestim


if perpeh==1{
matrix mao_mit_est=mao_mit
matrix sum_mit_est=sum_mao_mit
matrix modelmat_m_e=modelmat_m/nmiver 
scalar diestim_mi_e=diestim_mi
scalar nestim_mi_e=nestim_mi
}
else if perpeh==2{
matrix mao_mit_sl=mao_mit
matrix sum_mit_sl=sum_mao_mit
matrix modelmat_m_s=modelmat_m/nmiver 
scalar diestim_mi_s=diestim_mi
scalar nestim_mi_s=nestim_mi
}
else if perpeh==3{
matrix mao_mit_ot=mao_mit
matrix sum_mit_ot=sum_mao_mit
matrix modelmat_m_o=modelmat_m/nmiver 
scalar diestim_mi_o=diestim_mi
scalar nestim_mi_o=nestim_mi

}


}


*************************************************************************************
*	Final reports
*************************************************************************************

matrix mao_mit_tot=mao_mit_est+mao_mit_sl+mao_mit_ot
matrix list mao_mit_tot
matrix sum_mit_tot=sum_mit_est+sum_mit_sl+sum_mit_ot
*matrix list sum_mit_tot
matrix sum_mit_all=sum_mit_est\sum_mit_sl\sum_mit_ot\sum_mit_tot
matrix list sum_mit_all
matrix modelmat_m=modelmat_m_e,modelmat_m_s,modelmat_m_o
matrix list modelmat_m
matrix ones=J(1,7,1)
matrix modelmatt=ones*modelmat_m
matrix list modelmatt
di diestim_mi_e " " diestim_mi_s " " diestim_mi_o
di nestim_mi_e " " nestim_mi_s " " nestim_mi_o

*************************************************************************************
*	Totals
*************************************************************************************

/* Report Total */
matrix mao=sum_mit_all

matrix smao=mao[1...,3]
matrix bmao=mao[1...,2]
matrix nmao=mao[1...,1]

*matrix tmao=bmao

matrix cil_mao=J(4,1,0)
matrix ciu_mao=J(4,1,0)
forval i= 1/4{
matrix smao[`i',1]=sqrt(smao[`i',1])
*matrix tmao[`i',1]=bmao[`i',1]/smao[`i',1]
matrix cil_mao[`i',1]=int(bmao[`i',1]-signi*smao[`i',1]+0.5)
matrix ciu_mao[`i',1]=int(bmao[`i',1]+signi*smao[`i',1]+0.5)
}


matrix mao=nmao,cil_mao,bmao,ciu_mao,bmao,smao
matrix list mao

svmat mao, names(v)
ren v1 n
ren v2 cl
ren v3 e
ren v4 cu
ren v5 em
ren v6 s

export excel n cl e cu em s using "$excel/${name}smi.xls", sheet("totals") sheetmodify firstrow(varlabels) missing(".")
clear


/* Report Strata */
foreach k in mao_mit_est mao_mit_sl mao_mit_ot mao_mit_tot{
matrix mao=`k'

matrix smao=mao[1...,3]
matrix bmao=mao[1...,2]
matrix nmao=mao[1...,1]

*matrix tmao=bmao

matrix cil_mao=J(nstrata,1,0)
matrix ciu_mao=J(nstrata,1,0)
forval i= 1/58{
matrix smao[`i',1]=sqrt(smao[`i',1])
*matrix tmao[`i',1]=bmao[`i',1]/smao[`i',1]
matrix cil_mao[`i',1]=int(bmao[`i',1]-signi*smao[`i',1]+0.5)
matrix ciu_mao[`i',1]=int(bmao[`i',1]+signi*smao[`i',1]+0.5)
}


matrix mao=nmao,cil_mao,bmao,ciu_mao,bmao,smao
matrix list mao

svmat mao, names(v)
ren v1 n
ren v2 cl
ren v3 e
ren v4 cu
ren v5 em
ren v6 s

export excel n cl e  cu em s using "$excel/${name}smi.xls", sheet("s`k'") sheetmodify firstrow(varlabels) missing(".")
clear


}




/* Report regions */
run "$code/regions"
foreach k in mao_mit_est mao_mit_sl mao_mit_ot mao_mit_tot{
matrix mao=`k'

matrix mao=region'*mao
scalar nregion=rowsof(mao)

matrix smao=mao[1...,3]
matrix bmao=mao[1...,2]
matrix nmao=mao[1...,1]

*matrix tmao=bmao

matrix cil_mao=J(nregion,1,0)
matrix ciu_mao=J(nregion,1,0)
forval i= 1/8{
matrix smao[`i',1]=sqrt(smao[`i',1])
*matrix tmao[`i',1]=bmao[`i',1]/smao[`i',1]
matrix cil_mao[`i',1]=int(bmao[`i',1]-signi*smao[`i',1]+0.5)
matrix ciu_mao[`i',1]=int(bmao[`i',1]+signi*smao[`i',1]+0.5)
}


matrix mao=nmao,cil_mao,bmao,ciu_mao,bmao,smao
matrix list mao

svmat mao, names(v)
ren v1 n
ren v2 cl
ren v3 e
ren v4 cu
ren v5 em
ren v6 s

export excel n cl e cu em s using "$excel/${name}smi.xls", sheet("r`k'") sheetmodify firstrow(varlabels) missing(".")
clear


}

*************************************************************************************
*	Differences
*************************************************************************************

/* Report Total */
matrix mao=sum_mit_all
matrix dmao=mao[2,1..2]-mao[1,1..2]    /*d*/
matrix smao=sqrt(mao[2,3]+mao[1,3])   /*t*/
matrix tmao=dmao[1,2]/smao[1,1]    /*t*/
matrix pmao=1-normal(abs(tmao[1,1]))    /*p*/
matrix cil_mao=int(dmao[1,2]-signi*smao[1,1]+0.5)
matrix ciu_mao=int(dmao[1,2]+signi*smao[1,1]+0.5)
matrix smao=int(smao[1,1]+0.5)   /*t*/
matrix dmao=dmao,smao,tmao,pmao,cil_mao,ciu_mao
matrix list dmao

svmat dmao, names(v)
ren v1 n
ren v2 e
ren v3 s
ren v4 t
ren v5 p
ren v6 cl
ren v7 cu
format %7.4f  p
export excel n e s t p cl cu using "$excel/${name}smi.xls", sheet("dmao") sheetmodify firstrow(varlabels) missing(".")
clear


/* Report Strata */
matrix dmao=mao_mit_sl[1..58,1..2]-mao_mit_est[1..58,1..2]    /*d*/
matrix smao=mao_mit_sl[1..58,3]+mao_mit_est[1..58,3]   /*v*/

matrix tmao=J(nstrata,1,0)
matrix pmao=J(nstrata,1,0)
matrix cil_mao=J(nstrata,1,0)
matrix ciu_mao=J(nstrata,1,0)
forval i= 1/58{
matrix smao[`i',1]=sqrt(smao[`i',1])
matrix tmao[`i',1]=dmao[`i',2]/smao[`i',1]
matrix pmao[`i',1]=1-normal(abs(tmao[`i',1]))    /*p*/

matrix cil_mao[`i',1]=int(dmao[`i',2]-signi*smao[`i',1]+0.5)
matrix ciu_mao[`i',1]=int(dmao[`i',2]+signi*smao[`i',1]+0.5)
matrix smao[`i',1]=int(smao[`i',1]+0.5)   /*t*/
}

matrix dmao=dmao,smao,tmao,pmao,cil_mao,ciu_mao
matrix list dmao

svmat dmao, names(v)
ren v1 n
ren v2 e
ren v3 s
ren v4 t
ren v5 p
ren v6 cl
ren v7 cu
format %7.4f  p
export excel n e s t p cl cu using "$excel/${name}smi.xls", sheet("sdmao") sheetmodify firstrow(varlabels) missing(".")
clear




/* Report regions */
matrix mao_mit_sl=region'*mao_mit_sl
matrix mao_mit_est=region'*mao_mit_est

matrix dmao=mao_mit_sl[1..8,1..2]-mao_mit_est[1..8,1..2]    /*d*/
matrix smao=mao_mit_sl[1..8,3]+mao_mit_est[1..8,3]   /*v*/

matrix tmao=J(nregion,1,0)
matrix pmao=J(nregion,1,0)
matrix cil_mao=J(nregion,1,0)
matrix ciu_mao=J(nregion,1,0)
forval i= 1/8{
matrix smao[`i',1]=sqrt(smao[`i',1])
matrix tmao[`i',1]=dmao[`i',2]/smao[`i',1]
matrix pmao[`i',1]=1-normal(abs(tmao[`i',1]))    /*p*/

matrix cil_mao[`i',1]=int(dmao[`i',2]-signi*smao[`i',1]+0.5)
matrix ciu_mao[`i',1]=int(dmao[`i',2]+signi*smao[`i',1]+0.5)
matrix smao[`i',1]=int(smao[`i',1]+0.5)   /*t*/
}

matrix dmao=dmao,smao,tmao,pmao,cil_mao,ciu_mao
matrix list dmao

svmat dmao, names(v)
ren v1 n
ren v2 e
ren v3 s
ren v4 t
ren v5 p
ren v6 cl
ren v7 cu

format %7.4f  p
export excel n e s t p cl cu using "$excel/${name}smi.xls", sheet("rdmao") sheetmodify firstrow(varlabels) missing(".")

clear
capture log close
